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The dissipative Dicke model exhibits a fascinating out-of-equilibrium many-body phase transition 
as a function of a coupling between a driven photonic cavity and numerous two-level atoms. We 
study the effect of a time-dependent parametric modulation of this coupling, and discover a rich 
phase diagram as a function of the modulation strength. We find that in addition to the established 
normal and super-radiant phases, a new phase with pulsed superradiance which we term dynamical 
normal phase appears when the system is parametrically driven. Employing different methods, we 
characterize the different phases and the transitions between them. Specific heed is paid to the role 
of dissipation in determining the phase boundaries. Our analysis paves the road for the experimental 
study of dynamically stabilized phases of interacting light and matter. 
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Experimental progress in control and manipulation of 
light-matter quantum systems has generated a growing 
interest in many-body phenomena out of equilibrium [I] . 
Well established examples of such systems include ultra¬ 
cold atomic or ionic quantum gases in high finesse opti¬ 
cal cavities [2], semiconductor microcavities in the strong 
coupling regime mi, and superconducting qubits in mi¬ 
crowave resonators mm. The engineered interplay be¬ 
tween light and matter in these systems has led to the 
observation of a host of fascinating collective phases and 
quantum phase transitions including superfluidity in po- 
laritons [6] and the super-radiant Dicke phase transition 
in a Bose-Einstein condensate (BEC) coupled to an op¬ 
tical cavity [7] . 

A defining feature of many of these systems is that 
they are inherently driven and subject to dissipation. 
Driven dissipative systems are usually treated within a 
rotating frame formalism. This effectively renders the 
problem time-independent with the important feature 
that the asymptotic steady states are necessarily out of 
equilibrium. However, parametric driving of the system 
often does not allow the usual rotating frame simplifi¬ 
cations. Consequently, the resulting interplay between 
interactions, dissipation, and parametric driving, could 
lead to novel and exotic steady-state physics that has no 
counterpart in the undriven case. Parametric driving is 
increasingly used as an experimental tool in diverse con¬ 
texts, e.g. , in the generation of Floquet topological insu¬ 
lators 0, improved measurement fidelity with squeezed 
quantum states [3 HD] and unconventional phenomena in 
cavity QED [TT] . 

A prime example of a system exhibiting light-matter 
collective phenomenon is the Dicke model m • Here, a 
bosonic/cavity mode is coupled to a large number of two- 
level atoms. It exhibits a Z 2 symmetry breaking quantum 
phase transition from a normal phase (NP), where all the 
atoms are in their ground state and the cavity is empty, to 
a super-radiant phase (SP), where the atoms are excited 
and the cavity is in a coherent state. This model has 
recently been realized by coupling the external degree 


of freedom of a BEC to a quantized mode of a laser- 
driven optical cavity, and the theoretically predicted non¬ 
equilibrium phase transition has been observed 0 EES. 
Moreover, the inevitable photon leakage out of the cavity 
as well as dissipation of the BEC has been shown to lead 
to a considerable modification of the critical exponents 
of the transition DU- 



FIG. 1. A sketch of a parametrically driven Dicke model. A 
single-mode cavity is driven by a laser beam with an oscillat¬ 
ing laser-field power P(t). The cavity is naturally leaky with a 
dissipation rate n. Inside the cavity, an atomic cloud is cooled 
and forms a Bose-Einstein condensate that is coupled to the 
driving laser with coupling A (t) oc \/P(t) [cf. Eq. ([Tj)] . The 
atoms are also coupled to an environment with a dissipation 
rate r\. The system is best described by Liouvillian dynamics 
[cf. Eq. §]. 

In this work, we analyze the impact of parametric driv¬ 
ing on the phase diagram of the dissipative Dicke model. 
Specifically, we consider a modulation of the atom-cavity 
coupling, which is easily realizable in current experimen¬ 
tal setups, see Fig. [l] Using a combination of mean field 
theory and effective Hamiltonians, we obtain a rich phase 
diagram comprising: (i) the NP with parametric amplifi¬ 
cation, (ii) the SP phase, and (iii) a novel dynamical nor¬ 
mal phase (D-NP), which appears to be a dynamically ro¬ 
tating NP with pulsed super-radiance. We elucidate the 
vital role that dissipation plays in modifying the com¬ 
plex phase topography of this nonequilibrium system. 
Our analysis presents parametric driving as a promis¬ 
ing frontier in the search for exotic collective phases in 
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FIG. 2. Characteristics of the dynamical phase diagram of the parametrically driven Dicke model as a function of cavity/atoms 
frequency ujq — uj c — uo a and parametric modulation strength e for bare coupling Ao = 0.4Q and dissipation rates k — 77 = 0.1Q. 
We obtain these characteristics from the numerical steady-state solutions of the mean field equations [cf. Eqs. §-©]• Each 
point, with a 10 -3 resolution, is a result of a different numerical integration. Superimposed (striped overlay) is the normal 
mode stability zone [cf. Eq. §]. (a) Density plots of the absolute value of the steady-state time-averaged order parameters, (b) 
In steady-state the order parameters are oscillating with a complex beat structure m- The density plots are of the amplitude 
of the maximal frequency uj contributing to this oscillation. We see three distinct phases appearing as a function of e, i.e. , 
extensions of the normal and super-radiant phases (NP and SP), as well as a novel dynamical phase (D-NP). These three phases 
meet at a shared multicritical point. The properties of each phase are summarized in Table [I] 


light-matter systems. 

The single mode parametrically driven Dicke model is 
described by the Hamiltonian 

H(t) = tudca) a + fiiOg s l z + —j=^ yi 4 ( Q + at ), W 

i=l v iV 

where s l a with a = x,y, z are the spin operators describ¬ 
ing the i th two level atom, and a, a) represent the cavity 
creation and annihilation operators. The cavity’s reso¬ 
nance frequency is u; c , whereas the atoms are considered 
to be identical with level spacing hw a . We consider a 
time-dependent coupling between the atoms and the cav¬ 
ity of the form A (t) = Ao + e cos(2 fit). Such a coupling is 
easily generated by a modulation of the laser power that 
drives the cavity m ■ For e = 0 , the system exhibits the 
well known continuous phase transition from a NP to a 
SP when the coupling Ao > y/w c uj a /2 [ 12 ]. For e / 0 , 
we reiterate that the parametric driving described here 
cannot be rotated away by a suitable choice of frame. 
Indeed, recent treatments of similar modulations of the 
Dicke model were addressed using a mapping to para¬ 
metric oscillators, and a partial phase diagram for the 
NP was obtained pH EE]. Here, we explicitly include 
dissipation for both the cavity and the atoms (see Fig. [I]) 
and analyze the impact of parametric driving on the full 
phase diagram of the dissipative Dicke model, see Fig.[2j 

The driven and dissipative nature of the system is de¬ 
scribed by a Liouvillian equation for the density matrix 



Stability A\ 

Stability A2 

M 

\x\ 

\v\ 

1*1 

NP 

yes 

yes 

0 

0 

0 

1/2 

SP 

no 

yes 

^0 

^0 

^0 

^0 

D-NP 

no (★) 

no (★) 

0 

0 

0 

< 1/2 


TABLE I. Summary of normal mode [cf. Eq. ([ 5 |] and mean 
field [cf. Eqs. ([6])-([8])] analyses. In Fig. [ 2 ja) , we observe 
three main regions, dubbed normal phase (NP), super-radiant 
phase (SP), and dynamical normal phase (D-NP). Each region 
manifests a different behavior summarized here, i.e. , which 
normal mode is stable in each region m , and what are the 
values of the different order parameters. In the D-NP region, 
we denote by “no (★)” that at least one normal mode is un¬ 
stable. 

p sys of the system 

^ 5 Psys] T [2(2/) S y S (2^ {(2^(2,/? S y S }] 

N 

+ ] [2<sl/9 S y S <S^ — ,/9 sys }] , (2) 

i ,3 = 1 

where s± = s 3 x + is J y are ladder operators. The first 
term on the r.h.s. describes the standard Hamiltonian 
evolution and the last two terms represent the Marko¬ 
vian dissipation for both cavity and a global dissipation 
for the atoms in Lindblad form with rates k and 77, re¬ 
spectively PI- Note that this approach is valid in the 
Born-Markov limit of weak dissipation. 

In the absence of paramteric driving, e = 0 , the NP is 
well described by considering the collection of two-level 
atoms as constituting a giant spin S = JN s z aligned 
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along the z axis [20]. The deviations of this giant spin 
away from this quantization axis can be characterized by 
the standard Holstein Primakoff representation for the 
spin operators, S z = b^b = y/N — tfbb, and 

S+ = b^y/N — tfb, where 6 , b t are standard bosonic oper¬ 
ators [20 . This approach can be extended to address the 
stability of the NP in the presence of parametric driving, 
e^O. Since deviations of S from the z axis are expected 
to be small, as N —>> oo we can map the Dicke Hamil¬ 
tonian [cf. Eq. ([l])] onto the problem of two harmonic 
oscillators whose coupling is parametrically driven, 

i^NP (t) = huj c a^a + huj a b^b + HX(t)(b + b^)(a + at ). (3) 

Focusing on the case cj 0 = = cj c , Eq. <© can be 

diagonalized in terms of normal modes, 

H^p(t) = h£l\(t)A^A\ T Ml 2 it)A^A2 , (4) 

where for m = 1,2, fl^(t) m uj q — (—l) m 2A(t)u;o are 
time-dependent normal mode frequencies, and A m = 
[<5+(a - (-l) m 6) +5-(a t - (-l) m 6 t )] are the corre- 
sponding normal mode operators with coefficients = 
2 v^o° . For computational simplicity, we also assume 

Each normal mode is a quantum Mathieu paramet¬ 
ric oscillator, i.e. , its fundamental harmonic frequency 
varies sinusoidally in time [22j [23]. The stability of 
each quantum parametric oscillator can be deduced from 
its displacement Tr {p(t)(A m + A^)}. It results in a 
complex stability diagram comprising “Arnold tongues” 
which delineate regions where the displacement, though 
parametrically amplified, remains bounded (stable), and 
those where the displacement grows exponentially with 
time (unstable) [15] 24). 25]. Incidentally, the result¬ 
ing stability diagram for the quantum oscillator is the 
same as that of the classical damped Mathieu oscilla¬ 
tor obeying the classical equations of motion for the 
displacement |15j |24] [25] , 

Xm +'jx m + Vt 2 m (t)x m — o, (5) 

The combination of the stability diagrams of the two 
normal modes yields the stability of the NP, i.e. , the NP 
is stable only if both dissipative normal modes are stable, 
see shaded area in Fig. [2] In the absence of dissipation, 
k = r] = 7 = 0, each normal mode A rn is unstable in 
the limit of infinitesimal parametric driving e —>• 0 at res¬ 
onant frequencies / m?n = yj (Ao) 2 + n 2 Q 2 + (—l) m (Ao) 
where n = 1,2,3... [ 22 ] [23]. We find that the lower 
boundary of the NP is dictated by the lowest Arnold 
tongue of A 2 , i.e. , where the time-independent part of 
1^2 (£) becomes negative. The remaining stability bound¬ 
aries of NP are determined by frequencies where either 
A m becomes unstable. The impact of disspation on the 
stability of NP can be understood from the physics of 


Mathieu oscillators where dissipation results in a modifi¬ 
cation of the stability criterion for the parametric oscilla¬ 
tor. In particular, weak Markovian dissipation leads to a 
pronounced stabilization of the NP in the vicinity of these 
resonant frequencies for small driving e Ao, and barely 
affects the stability at higher drive amplitudes EH [25]. 
Indeed in Fig. [2] we see substantial stabilization of NP 
at the resonant frequency Such stabilization of the 
NP in the many-body context of the Dicke model is a 
manifestation of the explicitly dissipation dependent non- 
equilibrium asymptotic state. 

From Fig. [2] we see that the NP occupies only a small 
part of the phase diagram when the system is paramet¬ 
rically driven. However, what lies beyond these stable 
NP regions cannot be accessed by the current approach 
and requires another method, such as mean field theory, 
which is well justified for the Dicke model in the limit of 
N 1. This method was successfully used for studying 
the SP which has broken Z 2 symmetry, for the non-driven 
case e = o nanoi nu. The mean field ansatz that we 
use states that the total density matrix in the steady 
state is a product state of the individual density matri¬ 
ces, Psys = Pc ® nZi <8>Pi, where p c and pi are density 
matrices of the cavity and the i th atom, respectively m • 
Furthermore, since all atoms are identical, we assume all 
pi to be equivalent. Substituting this ansatz into Eq. 
we obtain a set of coupled non-linear equations for the 
mean-field order parameters of the parametrically driven 
and dissipative Dicke model 


a = — iuo c a — 2i\{t)x — Ka , 

(6) 

x = - io a y - rjx , 

(7) 

V =u a x ~ 2A(t) [a + a*]z-r)y, 

(8) 


where we have defined the order parameters a = 
(a)/VN, x = (T,iSf)/N, y = (Ei4)/V z = 
(Si S ?)/V an d have assumed z = \J (1/2) 2 — |a;| 2 — |y| 2 
since the mean field equations satisfy the constraint 
x 2 + y 2 + z 2 = 1/4. 

Equations form a set of coupled non- 

autonomous differential equations. Using numerical stiff 
ordinary differential equations solvers, we integrate this 
set of equations to a long time limit, where we ob¬ 
tain a convergent behavior. We find that, generically, 
the steady-state mean-field solutions show oscillatory be¬ 
havior around a zero or non-zero mean value [15]. In 
Figs- Up), we plot the absolute time-averaged order pa¬ 
rameters, |a|, |x|, \y\ and |z|, averaged over a sufficiently 
long time window in the steady-state[15]. Note, that the 
absolute value is taken for presentation reasons only, we 
always have a,y,z > 0 and x < 0. Based on these mean 
values, we find that the SP, characterized by a ^ 0, ex¬ 
tends from its zero drive region of ujq < 2Ao onto a large 
regime spanning both small and large drive amplitudes. 
At the critical frequency uj cr [ t (cjo , e), we observe a tran¬ 
sition to a region with a = x = y = 0. Contrasting these 








4 


results with those obtained through the study of normal 
modes, we see that for e < 0.3S4, the NP lies above the 
line defined by a; cr i t with z = 1/2. The NPe^SP tran¬ 
sition in this regime is thus an extension of the usual 
continuous Dicke transition at zero drive to finite para¬ 
metric driving. Note that the details of the transition 
may still differ from the standard Dicke transition as the 
parametrically-driven NP accommodates a large number 
of photons in the cavity. 

Interestingly, at e ~ 0.3S4, we see a sudden change 
in the curvature of cj C rit- This exactly signals the point 
where the NP ends and a novel dynamical phase, dubbed 
dynamical-NP (D-NP), starts. As opposed to the NP, 
though a = 0, this phase has oscillatory a(t) and does not 
have its spin aligned along the z-axis, i.e. ,z< 1/2. Addi¬ 
tionally, this region corresponds exactly to the parametri¬ 
cally unstable Arnold tongues of the aforementioned nor¬ 
mal modes. As a result, the point (e ~ 0.3f2, cj cr i t ), ap¬ 
pears to be a multicritical point where the three phases: 
NP, SP, and the new D-NP intersect. The principal fea¬ 
tures of the three phases are summarized in Table [T| 



FIG. 3. Trajectories on the Bloch sphere of the atoms 
order parameters as a function of time in steady-state, 
f[x(t),y(t),z(t),t\. The trajectories are for e = 0.50, Ao = 
0.40, and k = rj = 0.10. The trajectory that does not en¬ 
circle the z -axis (magenta) is in SP with ujq = 0.60. The 
trajectories that encircle the z -axis (red and yellow) are for 
ujq — 0.650 and ujq — 1.50, respectively. 

To better understand the nature of the D-NP, as well 
as the effects of the drive e / 0 on the NP and SP, 
we analyze the oscillatory behaviour around the steady- 
state mean-field solutions of Eqs. Typically, we 

find that each phase has a different oscillatory behav¬ 
ior: (i) in NP, the order parameters converge to zero and 
do not oscillate, (ii) in SP, the oscillations are small but 
mildly grow with e, and (iii) in D-NP, the order param¬ 
eters oscillate strongly around zero m We, then, Fast 
Fourier Transform (FFT) the steady-state solutions for 
each cjo and e. In both SP and D-NP regimes, the oscilla¬ 
tions have a complex beat structure that corresponds to a 
comb of peaked frequencies at uo ^ 0, as well as uo ^ £4, c^o 
m- Note that the frequencies that appear can be un¬ 


derstood from Floquet analysis [27.. In Figs. [2|b), we 
plot the amplitude of the largest uj ^ 0 peak in the FFT 
landscape in order to quantify the overall extent of the 
oscillation. We find that, SP has weak oscillations in all 
of the order parameters, whereas the D-NP is strongly 
oscillatory in the x — y plane. These aspects are better 
highlighted in Fig. [3j by plotting the steady-state time- 
dependent trajectory of the total spin (J2 i s' l \(t)) on th e 
Bloch sphere for different parameter values. It appears 
that a distinguishing criterion between the SP and D-NP 
is whether the trajectory encircles the z- axis. Our re¬ 
sults seem to indicate that the most plausible candidate 
for the D-NP is a ” normal phase” in a dynamical rotating 
frame. 

Combining the results from the normal modes and 
mean field analyses, we see that periodic modulation of 
the atom-light coupling results in a rich phase diagram, 
characterized by a multitude of dynamical phase bound¬ 
aries between the NP, SP and the intriguing new phase, 
D-NP. All three phases meet at a multicritical point. In 
the D-NP, the cavity periodically emits pulses of pho¬ 
tons with opposing phases, which should be detectable 
experimentally. The NP —>> SP boundary is principally 
dictated by where the normal mode A 2 becomes unsta¬ 
ble, whereas the NP —>■ D-NP boundary is fixed by the 
instability of any of the modes A m [15] . Within our mean 
field approach, we find the transitions SP —> NP, D-NP to 
be continuous, though the latter is rather sharp. How¬ 
ever, the nature of the NP —>• D-NP transition cannot 
be studied within our approach. The topography of the 
phase diagram is expected to vary with the choice of Ao 
and the strength of dissipation. Consequently, though 
the three phases would exist, the highly sensitive multi¬ 
critical point may disappear. 

Remarkably, we see that dissipation leads to a sizable 
stabilization of the NP. Due to the parametric nature of 
the normal modes, the NP in the driven case can manifest 
a dissipation-assisted generation of substantial entangle¬ 
ment/squeezing between the atoms of the condensate and 
the cavity [21]. The physical signature of such entangle¬ 
ment as well as the impact of parametric driving and 
dissipation on the critical exponents defining the differ¬ 
ent phase transitions merit in-depth studies. It would 
also be interesting to extend the present work to other 
parameter regimes like uj a « uo c realized in current ex¬ 
perimental setups [7j. 

Our work shows that parametric driving is a powerful 
tool in the quest for new physics, which exists exclu¬ 
sively in the realm far from equilibrium. The richness 
of the physics seen in the simple Dicke model presages 
intriguing phenomena in time-dependent systems, which 
requires the development of new theoretical methodolo¬ 
gies. This frontier is potentially best explored using ex¬ 
perimental light-matter systems. 
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SUPPLEMENTAL MATERIAL 


I. NORMAL PHASE 


A normal mode with time-dependent frequency mod¬ 
ulation [cf. Eq. @ in the main text] can be generically 
described by the Hamiltonian for a parametric oscillator: 

H =^n + \H+n{t)W- (LI) 

For a classical oscillator with sinusoidal modulation, the 
above Hamiltonian describes the well known Mathieu 
parametric oscillator [22l EH- The quantum parametric 
Hamiltonian can be rewritten as 

H = [cj 0 T ^-^-]o)a + ^-^(a 2 + a t2 ) 

2i(jUq 4co>o 

= \fFFF)a ] a, (1.2) 

where the ladder operators a, at are defined with respect 
to the time independent model (i.e., fi(t) = 0). 



FIG. 1.1. Characteristic numerical time-integration plots of 
Eq. (1.8) in different regions of the phase diagram. All curves 
are with Ao — 0.49, e — 0.59, and k = r\ — 0.19. We see 
that for (i) uo = 0.50, mode A\ is stable whereas A 2 is not, 

(ii) coo = 0.70, mode A\ is unstable whereas A 2 is stable, and 

(iii) coo = 0.90, both modes A m are stable. 


We consider a coupling to an external bath which is in 
the rotating wave approximation. Hence, the solutions 
to the Heisenberg equations of motion, in the presence of 
dissipation, for the operators a and aJ take the form [28] 

a(t) = G(t)a(0 ) + L* (; t)a 1 (0) + F{t ), (1.3) 

o) (t) = G*(t)a^( 0) + L{t)a{0 ) + F^(t ), (1.4) 

where G and L are time dependent functions obeying the 
initial conditions G(0) = G*(0) = 1 and L(0 ) = L*(0) = 
0. F is an operator term which stems from the dissipation 
and also depends on the functions G and L. It satisfies 
the condition F( 0) = U(0) = 0. The functions G, L 
obey the integro-differential equations 

G(t) = -i(wo + ^)G - i^L - f dsK{t - s)G(s) , (1.5) 

ZidQ /(do Jo 

L(t) = i(uoo + — [ dsK(t — s)L(s ), ( 1 . 6 ) 

2ojo 2ojq Jq 


Tr{p(q ( a x + 4) } Tr{p(d ( a 2 + } 
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FIG. 1.2. The numerical stability diagram of the normal 
modes A m [cf. Eq. in the main text] for Ao = 0.49 and 
k = rj = 0.19. As the drive in the Dicke model affects each 
mode differently, we obtain two “Arnold tongue” stability di¬ 
agrams that are shifted and scaled with respect to each other. 
Superposing the zones where both modes are stable leads to 
the stability of the normal phase (NP) [cf. Fig. [5] in the main 
text]. 


where the dissipative kernel K(t) = f duuJ(uj)e~' luJt and 
J(uj) is the spectral density of the dissipative bath. 

These coupled first order integro-differential equations 
are rather difficult to solve and numerical solutions are 
needed. However, for standard Markovian dissipation, 
induced by cavity leakage in the rotating frame or cou¬ 
pling to an ohmic bath, K(t ) = 7 S(t) where 7 is the 
damping rate. Substituting this in Eqs. (1.5) and (1.6), 


we see that they become a set of linear ODEs. Observ¬ 
ables and correlation functions can easily be obtained 
from these solutions. For example, 


<«(*)> = G(t)(a( 0)> + L*(t)(a\0)) + (F(t)), (1.7) 


where the expectation values are with respect to the ini¬ 
tial density matrix. Assuming initial conditions such that 
(. F(t )) = 0, which are expected for baths with no partic¬ 
ular ordering, we obtain 

= Tr {P(t) (a + a f )} (1-8) 

= R e[G(t) + L(t)]{x{ 0)) - lm[G(t) + . 

muj o 

For arbitrary initial conditions, the stability of the oscil¬ 
lator is dictated by whether the pre-factors [G(t) + L(t)\ 
grow with time as one approaches the asymptotic state. 
For the parametric oscillator, we expect it to become 
exponentially unstable as the strength of the driving is 
increased m- Choosing fi(t) = gcos(29t), the zones of 
stability can be traced in the ujq — g plane. The result¬ 
ing stability diagram is the same as that for the classical 
Mathieu oscillators, which can also be extracted from 
the classical equations of motions [cf. Eq. © in the main 
text]. 

To obtain the full NP stability diagram of the Dicke 
model [see Fig [2] in the main text], we study the stability 
of both normal modes A m , m = 1,2, with frequencies 
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FIG. II. 1. Characteristic numerical analysis of the mean-field 
equations [cf. Eqs. in the main text] leading to the dis¬ 

played phase diagram [cf. Fig. [ 2 ] in the main text]. All plots 
are with Ao — 0.4Q, e — 0.5S2, and k = 77 = O.lfT (a) and 
(c) are characteristic numerical time-integration plots of the 
mean-field equations with c^o = 0.6Q and ujo = 0.65Q, respec¬ 
tively. (b) and (d) are the corresponding Fast Fourier Trans¬ 
form (FFT). We see that in the super-radiant phase (SP) 
region [Figs, (a) and (b)], all order parameters are oscillat¬ 
ing around a non-zero mean with relatively small oscillations 
amplitudes. In the dynamical normal phase (D-NP) region 
[Figs, (c) and (d)], apart from z, all order parameters have 
a zero mean value, but their oscillations are large taking the 
full length of the central spin in its x and y components. 


[cf. Eq. Q in the main text] 

f2^(£) = cUq T 2AoCJo T 2iUJq€ cos(2f2t), (1.9) 

£l\(t) = cjq — 2AoCo>o — 2cjoecos(2f2t). (1.10) 

To simplify our calculation, we also assume that the 
baths the two modes couple to, have identical spectral 
densities [21]. Though relaxing this condition would lead 
to more technical complexity, it should not have any non¬ 
trivial physical consequence in the limit of weak dissipa¬ 


tion studied here. 

We solve the corresponding Eqs. (1.5)-(1.6) and in 
Fig. |I.1[ we present characteristic numerical time- 
integration plots of Eq. (|I.8|) for the normal modes A m . 


Repeating this procedure for different ojq and e, and by 
checking for converging/diverging solutions [cf. Fig 
we find the stability diagram for both modes A m 


1.1 


see 


Fig. |I.2| The superposition of the two stability diagrams 
then yields the stability of the normal phase shown in 
Fig. [2] 


II. MEAN-FIELD ANALYSIS 


In the main text, we obtained a set of coupled non- 
autonomous mean-field equations [cf. Eqs. in the 

main text] for the mean field parameters a,x and y. 
These equations were solved numerically with a variety of 
ODE solvers. Characteristic numerical time-integration 
plots of the solutions to these equations are shown in 
Figs. II.l| (a) and (b). Note that the solutions converge to 
the asymptotic regime for times tQ ~ 500 and are oscilla¬ 
tory. The time scale for reaching the asymptotic regime 
varies with the parameters. Repeating this procedure 
as a function of ujq and e, the time-average of the order 
parameters over the steady-state behaviour (the last one- 
third of the integrated time) is presented in Figs.^a) in 
the main text. 

We find that typically the solutions show sinusoidal 
oscillations characterized by the frequency of the para¬ 
metric drive. However, in certain parameter regimes, 
the order parameters oscillate strongly with a complex 
beat structure involving multiple frequencies. To ana¬ 
lyze all these solutions in a systematic manner, we Fast 
Fourier Transform (FFT) the steady-state signals, see 
Figs. [IT3 c) and (d). The largest amplitude of a finite 
frequency in such plots serves as a measure for the ex¬ 
tent of the oscillation that the order paramaters undergo. 
This amplitude is plotted as a function of ujq and e in 
Figs. §b) in the main text. 



































